#!/usr/bin/octave 

sizes = [35, 50, 75, 100, 139, 150, 160, 200, 250, 300, 340, 350, 400, 450, 490, 500];

load standard.txt
load sample.txt

height = max(standard(floor(end/3):end)); # discards leading garbage by picking 
                                          # the highest peak in the last two-thirds of the standard
quantum = height/10
threshold = 2*quantum*abs((sort(floor(standard ./ quantum)))(floor(end/2)));
if threshold == 0
  threshold = 100;
endif
threshold

x = 1:length(standard);
peak = find(standard > threshold);
first = peak(1) + 100; # the first peak is usually crap
lastpeak = peak(end);
last = lastpeak + 100; # don't want to mangle the last peak
[first, last]
% trim the leading garbage and negative values
outliers = find(standard > 1.5*height);
standard(outliers) = height;
outliers = find(standard < -threshold);
standard(outliers) = 0;

outliers = find(sample > 1.5*height);
sample(outliers) = height;
outliers = find(sample < -threshold);
sample(outliers) = 0;

plot(x(first:last), standard(first:last), x(first:last), sample(first:last));
pause;
#exit

lastpeak = lastpeak - first;
scale = max(sizes)/lastpeak;
size(standard)
standard = standard(first:last);
sample = sample(first:last);
x = (1:length(standard)) * scale;
standard(1:100) = 0; # zero out the crud
sample(1:100) = 0;

ref = zeros(1, length(standard));
for i = (1 : length(sizes))
  ref = ref + height * pulse ( sizes(i), 0.2, x ); 
endfor

#axis([110,370,-200,height]);
#plot(x, standard, "*", x, sample, "@", x, ref);
#plot(x, standard, "*", x, ref);
#pause;
#axis;

order = 3;
dwt = repmat(Inf, 1, 31);
for window = 5:2:31
  sm = sgolayfilt(sample, order, window);
  dwt(window) = dw(sample, sm);
  if dwt(window) < 1
    break;
  endif
endfor

[d, window] = min(abs(dwt - 2));
sm = sgolayfilt(sample, order, window); # redo smoothing with the optimal window

order = 2; window = 7;
d1 = sgolayfilt(diff(sm), order, window);
d2 = sgolayfilt(diff(d1), order, window);
d3 = sgolayfilt(diff(d2), order, window);

x = x(3:length(d3)+2);
sm = sm(3:length(d3)+2);
d1 = d1(2:length(d3)+1);
d2 = d2(2:length(d3)+1);

noise0 = 10*noise(sample);
noise1 = 10*noise(d1);
noise2 = 10*noise(d2);


#axis([250,300,-200,height]);
axis([14,19]);
#axis([100,200]);
#plot(x, sm, "*", x, d1, "g*-", x, d2, "b*-", x, d3, "o*-");
#plot( x, d2, x, noise2*ones(1,length(x)), x, -noise2*ones(1,length(x)));
#plot( x, sm, x, median(standard) + noise0*ones(1,length(x)), x, median(standard) - noise0*ones(1,length(x)));
#plot( x, sm, x, 10*d2, x, -10*noise2*ones(1,length(x)));

peaks = zeros(1, length(x));
region = find(sm > noise0 & (d1 > noise1 | d1 < -noise1 | d2 < -noise2));
peaks(region) = ones(1, length(region));

starts = find(diff(peaks) > 0);
ends = find(diff(peaks) < 0);
length(starts)
length(ends)

height = zeros(1,length(starts));
area = zeros(1,length(starts));
pos = zeros(1,length(starts));
for i = 1:length(starts)
  from = starts(i);
  to = ends(i);
  [height(i), p] = max(sm(from:to));
  pos(i) = x(from + p);
  area(i) = sum(sm(from:to));
endfor

format short g;

out = horzcat(pos', height', area')
file = fopen("out", "w");
for i = 1:length(out)
  fprintf(file, "%g\t%g\t%g\n", pos(i), height(i), area(i));
endfor
fclose(file);

for b = 0:20:480
  axis([b, b+20]);
  plot( x, sm, x, 5*d1, x, 4*d2, x, 5*d3, x, -5*noise1*ones(1,length(x)), x, 5*noise1*ones(1,length(x)), x(region), sm(region), "*", x, [50*diff(peaks),0], "b@");
  pause;
endfor
exit;

from = find(x > 110)(1)
to = find(x > 370)(1)
X = vertcat(ref, standard', sample');
X = X(:, from:to);
save X.txt X

